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ABSTRACT 

We have combined the thermo-chemical disc code ProDiMo with the Monte Carlo radiative 
transfer code MCFOST to calculate a grid of ^300000 circumstellar disc models, systemat- 
ically varying 1 1 stellar, disc and dust parameters including the total disc mass, several disc 
shape parameters and the dust-to-gas ratio. For each model, dust continuum and line radiative 
transfer calculations are carried out for 29 far IR, sub-mm and mm lines of [OI], [CII], ^^CO 
and 0/P-H2O under 5 inclinations. The grid allows to study the influence of the input param- 
eters on the observables, to make statistical predictions for different types of circumstellar 
discs, and to find systematic trends and correlations between the parameters, the continuum 
fluxes, and the line fluxes. The model grid, comprising the calculated disc temperature and 
chemical structures, the computed SEDs, line fluxes and profiles, will be used in particular for 
the data interpretation of the Herschel open time key programme GASPS. The calculated 
line fluxes show a strong dependence on the assumed UV excess of the central star, and on the 
disc flaring. The fraction of models predicting [OI] and [CII] fine-structure lines fluxes above 
Herschel/Pacs and Spica/Safari detection limits are calculated as function of disc mass. 
The possibility of deriving the disc gas mass from line observations is discussed. 

Key words: astrochemistry; stars: formation; circumstellar matter; radiative transfer; meth- 
ods: numerical 



1 INTRODUCTION 

The structure, composition and evolution of protoplanetary discs 
are important comer-stones to unravel the mystery of life, as they 
set the initial conditions for planet formation. Spectral energy dis- 
tributions (SEDs), although their analysis is known to be degen- 
erate, probe the amount, temperature and overall g eometry of the 
dust in the discs, such as disc flaring ('Meeus et al." 2001), puffed- 
up inner rims (DuUemond et al. 2001; Acke et al. 2009), and indi- 
cations of an average g rain growth in discs as young as a few Myr 
jP' Alessio et alj2001l) . 

Most works in the past dec ade have focused on the analysis of 
SEDs of individual objects (e.g. lD' Alessio et al]|200 (^ or to study 
systematic trends in infrared colours of various types of discs, rang- 
ing from embed ded young stellar objects to exposed TTauri stars 
( iRobitaille et alii2006i ). Some ambiguities inherent in SED analy- 
sis ca n be resolved by images in scatte red light ((sta pelfeldt et al. 
Il998h . in mid-infrared thermal emission jMcCabe et aU2003l) or in 



the mm-regime jAndrews & Williamsll2007l) . The Spitzer observa- 
tory has enabled detailed studies on dust mineralogy, constraining 
dust properties in the up per layers of the inner disk regions by us - 
ing solid-state features iFurlan et al. I I2OO6I ; IOlofsson etai]|2009l) . 
Multi-technique, panchromatic approaches, combining the afore- 
mentioned observations, are now becoming possible but remain 
limited to a few objects with complete dat a sets (e.g. iWolf et al.l 
l2003l ; |Pinte et al.ll2008l ; rDuchene et alj2009l) . 

However, all these observational findings are related to the 
dust component. Initially, about 99% of the disc mass can be as- 
sumed to be present in form of gas, and the progress toward a bet- 
ter understanding of the gas component, such as chemical com- 
position, gas temperature structure and vertical disc extension, is 
hampered by a current lack of observational data. With several 
key programs of the HERSCHEL SPACE OBSERVATORY, such as 
GASPS and WISH, a new body of gas emission line observa- 
tions will be provided very soon, and there is a clear need to 
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Figure 1. Unified MCFOST / ProDiMo model pipeline to calculate one radiation thermo-chemical disc model with SED and line flux prediction. 



include the gas in such systematic studies. Recent work has fo- 
cused on the prediction of far IR line emissions from individ- 
ual discs (Meiierink et alJl20()8l : rErcolano et Zll2008l : lwoitke et al.l 



Table 1. Parameters of the DENT grid and values assumed. Rsubu stands 
for the dust sublimation radius (where = 1500 K). The choice of incli- 
nation angles resembles a randomly oriented sample. 



2009; Cernicharo et al. 



20091). or rather small par ameter studies 



( Kamp et alJ l2009h . lG'oicoechea & Swinvarj ( |2009|) have recently 
discussed the detection rates of the [OI] 63 ^m, [SI] 56 fim and 
[Sill] 34 fine-structure lines with Herschel and the proposed 
SpiCA/Safari mission on the basis of one low mass disc model 
(Mgas = 10"'' M0) by[Gorti & Hollenbach (2004). Their conclu- 
sions, however, depend on the choices of the other model parame- 
ters, and it is difficult to put them on a firm statistical basis. 

To address these issues, we have combined our state-of-the- 
art computer codes to calcul ate the dus t and line radiative trans- 
fer with the MCFOST-code jPinte et al. 1 12006), and the gas ther- 
mal balance and chemis try with the PRODlMo-code ( I Woitke et al.l 
l2009l ; lKamp et alj2009l) . We have computed a large grid of 300 000 
disc models to simultaneously predict SEDs and gas emission lines 
from parametrised disc density distributions. The grid name DENT 
stands for Disc Evolution with Neat Theory. The following sec- 
tions describe the model pipeline (Sect. |2ll and the first results on 
fine structure emission line fluxes and detectability (Sect.[3}- Sec- 
tion |4] summarises the preliminary findings of the DENT grid and 
provides an outlook to future studies. 



2 THE DENT GRID 

The DENT grid is a tool to investigate the influence of stellar, disc 
and dust properties (Table [TJ on continuum and line observations 
(Table |2l(, and to study in how far these dependencies can be in- 
verted. The grid is designed to coarsely sample the parameter space 
associated with young, intermediate to low-mass stars at ages (1- 
30Myr) having gas-rich and gas-poor discs. 

In the DENT grid, 11 variable and 2 fixed input parameters 
(see Table [TJ are required to specify the star -1- disc systems. The 
input stellar parameters are mass and age. The corresponding 
effective t emperature s and l uminosities are from the evolutionary 
models of lSiess et alj ilOOd) . For the photospheric spectra, we use 
Kurucz stellar atmosphere models of solar abundance and match- 
ing Teff and \og{g). Because some young stars show significant 
accretion and/or chromospheric activity, an extra UV component is 
added to the spectrum in DENT. This extra UV component has an 
important impact on the disc chemistry and temperature. It is de- 
fined as /uv = iuv /Li, where Luv = Jgi^^^ Lx dX is the UV 
luminosity with assumed spectral shape L\ oc A"'^. 

There are other sources of energy in the discs. The interstel- 
lar irradiation is a ssumed to be isotro pic and fixed throughout the 
grid by (Eq. 27 of lWoitke et ai]|2009l, with x'^'^= !)• The cosmic 
ray ionisation rate of H2 is set to ^cr ~ l.TxlO"^^ s~^. We further 



stellar parameter 


AU stellar mass [Mg] 




0.5, 1.0, 1.5,2.0, 2.5 


age age [Myr] 




1,3, 10, 100 


/uv excess UV /uv = Luv/L* 




0.001,0.1 


disc parameter 


Afj disc dust mass [Mq] 


10-^ 


lO"'', 10-5, 10"", 10-3 


Pd/Ps dust/gas mass ratio 




0.001,0.01,0.1, 1, 10 


/Jin inner disc radius [-Rgubli] 




1, 10, 100 


Rout outer disc radius [AU] 




100, 300, 500 


e column density Nn{r) oc r~ 


-e 


0.5, 1.0, 1.5 


Ho scale height [AU] 


f 


fixed: 10 @ ro = lOOAU 


P discflai-ingi/(r) = Ho(;^ 


0.8, 1.0, 1.2 


dust parameter 


"mill minimum grain size [pm] 




0.05, 1 


Qmax maximum grain size [pm] 


s/2 


fixed: 1000 


s settling H{r, a) cx H{r) a~ 


0, 0.5 


inclination angle 


I 0° (face-on), 41.4°. 60°, 75.5°, 90 


' (edge-on) 



Table 2. List of output quantities: monochromatic continuum fluxes -Fcont 
and integrated, continuum-subtracted line fluxes -Fiinc- 



observable wavelengths [pm] 



Kiont 57 A-points between 0.1pm and 3500/im 

01 63.18, 145.53 

CII 157.74 

i2cO 2600.76, 1300.40, 866.96, 650.25, 520.23, 433.56, 371.65, 

325.23, 289.12, 260.24, 144.78, 90.16,79.36, 72.84 

0-H2O 538.29, 179.53, 108.07, 180.49, 174.63, 78.74 

P-H2O 269.27, 303.46, 100.98, 138.53, 89.99, 144.52 



assume that all discs are passive, i.e. there is no viscous heating. 
This will have an impact on the inner structure and temperature of 
some discs. However, the SED of several sources in t he literature 
are w ell fitted without viscou s heating, e.g. IM L upi jPinte et al] 
I2OO8I) and IRAS 04158+2805 dOlauser et al1l2008h . 

The density structure of the disc is parametric, with power- 
laws used for the surface density, the scale height, and the 
dust size distributio n. Optical properties of astronomical silicate 
( iDraine&Le^l 19841) are used to calculate the dust opacities. 

The chemical abundances are calculated selecting 9 atomic 
elements, 66 molecular and 5 ic e species, 950 reac tions, and ele- 
ment abundances as outlined in dWoitke et al]|2009h . The numeri- 
cal code used here feature an improved coupling between the UV 
radiative transfer and the calculation of the UV photo-rates (see 
iKamp et al.ll200^ . This version of DENT does not contain PAHs. 
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This will have an impact in particular on the Herbig AeBe's, with 
a possible increase of the gas temperature, especially in the upper 
layers of the discs. We do not expect a significant effect on the con- 
tinuum. The role of PAHs will be explored in a forthcoming paper. 

To calculate the DENT grid, two numerical codes were used in 
a sequence that we now describe (see Fig.[TJ. In phase 1, MCFOST 
solves the dust radiative transfer problem to obtain the dust temper- 
ature structure Td{r, z) and the internal mean intensities Jv{t, z). 
In phase 2, these data are transferred to ProDiMo which cal- 
culates the gas temperature structure Tg{r, z) assuming gas ther- 
mal balance, the chemical composition nmoi{r, z) assuming kinetic 
chemical equilibrium, and the level population of the gas species 
in the disc. An escape probability method is used to calculate the 
level populations. Phase 2 requires an additional gas parameter: the 
dust-to-gas ratio, pd/pg. Among the results of phase 2 are the to- 
tal species masses A/moi and averaged gas temperatures (r™°') for 
molG {O, C+, CO, H2O}, defined as 

,^mou _ /wn,oi(r,z)rg(r,z) dV 
J nnioi(r, z) dV 

where nmoi{r,z) is the particle density at position (r, z) in the 
disc. In phase 3, the level populations are transferred back to Mc- 
FOST to calculate the emission line profiles. The formal line trans- 
fer solutions are computed in 301 velocity bins on 100 x 72 par- 
allel rays organised in log-equidistant concentric rings in the im- 
age plane. A Keplerian rotation velocity field is assumed for the 
bulk disk kinematics, and a thermal + turbulent broadening with 
'^turb = 0.15 km/s is added. The calculations are completed by 
phase 4 running formal solutions of the dust continuum radiative 
transfer problem on the same rays. The calculated continuum and 
line intensities are post-processed to get the integrated line fluxes 
after continuum subtraction. Further details are listed in Table[2] 

In total, the DENT grid comprises 323020 disc models and 
SED calculations. A total number of 1610150 line flux calculations 
have been carried out for 29 spectral lines under 5 inclinations. 
We note that some parameter combinations may lead to unrealis- 
tic models, but they have been kept for the sake of completeness. 



3 RESULTS 

Figure |2] depicts all calculated line fluxes of [01] 63.2 /im in form 
of 3 histograms, underlining the depending on the assumed stel- 
lar UV excess. Models with high /uv (blue in the figure) have 
warm disc surfaces heated by the stellar UV with Tg ^ Td, and 
hence strong emission lines. Models with low /uv (red) have gas 
temperatures more equal to the dust temperatures, and hence less 
strong emission lines. The dependence of -Flinc on /uv is signifi- 
cant for all stars (also in other lines), but is less pronounced for Her- 
big Ae/Be stars which already produce photospheric soft UV, even 
for /uv = 0. For example, a Kurucz stellar atmosphere model with 
Teff = 8500 K, log(g) = 4 attains Luv ~ 0.097L* whereas a model 
with Teff = 4500 K, log(g) = 4 only attains Luv ~ 2 • 10"''L*. Ta- 
ble [3] lists the fraction of models predicting line fluxes of [01] and 
[CII] larger than (0.5h,3(T) detection limits of Herschel/Pacs 
and Spica/Safari at I40pc as function of disc gas mass. Note 
that, for massive discs, the line emitting regions may be optically 
thick in the continuum and may have Tg ^Td. In such cases, there 
is only very limited contrast between line and continuum, leaving a 
fair fraction of massive discs with non-detectable lines. 

Another clear trend in the models is depicted in Fig. [3] which 
shows the calculated line fluxes of [01] 63.2 /^m for a sub-selection 
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Figure 2. Dependence of line flux [OI] 63.2/im on the stellar UV excess 
/uv- The black histogram counts the DENT models that result in certain 
[OI] 63.2/^m fluxes at distance 140pc in 40 log-equidistant bins. The red 
histogram represents the low /uv = 0.001 models, and the blue dotted his- 
togram the high /uv = 0.1 models. The difference between high and low 
UV excess causes a difference of about 1 — 1.5 orders of magnitude in line 
flux. The arrows show the (3(T, 0.5h) detection limits of Spica/Safari 
and Herschel/Pacs, see dependence on disc mass in Table[3] 

Table 3. Fraction of DENT models predicting line fluxes larger than 
(0.5h,3(T) detection limits of Herschel/Pacs and Spica/Safari at 
140pc as function of total disc gas mass Mgas [Mq] and stellar UV ex- 
cess /uv- Each entity in this table is based on more than 5000 disc models. 
We neglect detection problems due to background confusion here. 



fv^''^"' 10"^ 10"*^ 10"^ 10-"* 10--' 10-2 


10-1 


l^[OI]63.2Mml > 1.2 X 10-17 W/m2 (HERSCHEL/PACS) 

0.1 0% 20% 51% 66% 70% 71% 71% 
0.001 0% 6% 17% 23% 26% 30% 34% 


l-f'loiiws.SMml > 4x10-18 w/m2 (Herschel/Pacs) 
0.1 0% 0% 9% 30% 45% 52% 
0.001 0% 0% 3% 7% 12% 18% 


57% 
23% 


l-f^[cnii57.7A.ml > 4x10-18 w/m2 (Herschel/Pacs 
0.1 0% 0% 17% 52% 56% 57% 
0.001 0% 0% 6% 14% 14% 14% 


56% 
13% 


l^[OI163.2Mml > 1.8x10-19 W/m2 (SPICA/S AFARI) 

0.1 65% 90% 96% 95% 96% 97% 
0.001 38% 65% 79% 83% 86% 86% 


96% 
85% 


l^^[OI1145.5,.ml > 1.2x10-19 W/m2 (SPICA/S AFARI) 

0.1 5% 44% 68% 80% 85% 87% 
0.001 1% 15% 29% 46% 56% 62% 


88% 
63% 


l^[CII1157.7Mml > 1.2x10-19 W/m2 (SPICA/S AFARI) 

0.1 0% 83% 95% 94% 95% 94% 
0.001 0% 53% 81% 78% 76% 73% 


93% 
69% 



of T Tauri disc models with high /uv as function of disc mass. At 
low disc masses, the discs are optically thin and the flaring of the 
disc (see definition of /3 in Table [T]! has only little influence on the 
line fluxes. However, with increasing disc mass, the inner disc be- 
comes optically thick and the computed line fluxes split up into two 
branches. For strong flaring (/? — 1.2), the fluxes of the emission 
lines steadily increase further, whereas for non-flaring discs (/3 ^ 1) 
they saturate at Afgas ^ 10"^ - lO—* Mq, henceforth called the 
"saturation disc mass". In other words, for massive T Tauri discs, 
high far IR line fluxes (e.g. [01] 63.2 ^m>3 x lO'i^ W/m^) are 
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Figure 3. Dependence of line flux of [OI] 63.2^m at distance d = 140 pc 
on the flaring parameter (3, as function of total disc gas mass Mgas . A sub- 
selection of 3456 TTauri models is plotted (Af* 1 Mq, age 1 Myr, 
/uv = 0.1, iJin = ^subli: s = 0, inclination angle ^ 60°, where the 
statistical range of line flux predictions due to the variation of the other 
input parameters in expressed by mean values and standard deviations. 

a safe indicator of disc (gas) flaring. In flared geometry, the disc 
surface is directly heated by the star, hence higher temperatures and 
stronger emission lines. In non-flared geometry, the surface layers 
are situated in the shadow casted by the dust in the inner disc re- 
gions, and the gas temperatures are cooler. In that case, a further 
increase in disc mass does not lead to stronger emission lines, but 
rather to an increase of the shielding effects, causing cooler tem- 
peratures and sometimes even weaker emission lines. 

The saturation behaviour of the emission lines for non-flared 
geometry depends on the line properties. High-excitation lines like 
[OI] 145.5 react more sensitively to temperature changes and 
hence to disc flaring, whereas low excitation lines like CO J=l— >-0, 
[CII] 157.7 /im are less affected. However, the point where these 
saturation effects start to appear, the saturation disc mass, is found 
to be rather similar for all lines, because it is the amount of dust 
and its opacity in the inner disc regions that causes the shielding. 

The fact that different emission lines originate in different disc 
regions, and the strong dependencies of the line fluxes on UV ex- 
cess /uv and flaring index /3 suggest that a simple, PDR-like anal- 
ysis of emission line ratios to determine the total disc mass is diffi- 
cult. However, Fig.|4]shows that most of these complicated param- 
eter dependencies affect the line fluxes in an indirect way, namely 
by changing the mean temperature of the disc. If we plot the depen- 
dency between disc gas mass and [OI]63.2/^m line flux for models 
with similar mean disc temperatures, we roughly retrieve a linear 
relation as expected from a simple analysis (see Appendix lAt. 

The Icr errorbars in Fig. |4] show, however, that there is still a 
considerable variation of the [OI]63.2/im flux among the models, 
even if they have similar mean disc temperatures. Note also, that 
the fitting value of Tcxc = 30. 5K is not consistent with the actual 
disc temperatures as measured from the models. This may be due to 
the way we have defined the mea n disc temperature ( Eq.[T) and/or 
due to included non-LTE effects. iKamp et alj hoOSl) have shown 
that, in general, Tcxc < Tg in the outer and upper disc layers. 

We note that an unmindful usage of Eq. ( IA3I ) with an assumed 
value for Tcxc (e.g. 27K) for the purpose of gas mass determi- 
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Figure 4. Relation between total disc gas mass and [OI] 63.2/im line flux 
at distance 140pc. The I'-axis is subdivided into 30 log-equidistant bins. 
In each bin, the values of the disc gas mass (parameter, see Table [T) of 
all models with matching [OI] 63.2/xm line flux are statistically analysed to 
yield a mean value and a standard variation, which are then plotted as points 
with veilical error bars. We distinguish between cool and warm models by 
(TJ-*' ) 2S defined in Eq. {!), and have chosen temperature interval bound- 
aries to bracket the Icr distribution as shown by the inserted histogram. The 
following models have been selected: dust/gas ratio Pd/Pg = 0.01 — 1, 
flaring index /3 ^ 1, and inclination angle ^ 60°, altogether 122269 cool 
and 101244 warm models. The full and dashed lines show the results from 
ttie formula Eq. jA3) with Tcxc = 24K and Tcxc = 30. 5K, respectively. 

nation from a measured [OI]63.2/im line flux can be misleading 
and does not account for the variety of results that we find in the 
DENT model grid. In particular, high mass discs tend to be cooler 
as demonstrated by the saturation behaviour depicted in Fig. [3] and 
their oxygen fine-structure lines can easily become optically thick. 



4 SUMMARY AND CONCLUSIONS 

In a concerted effort of the theory groups in Edinburgh, Grenoble 
and Groningen, we have computed a grid of 300 000 circumstellar 
disc models, simultaneously solving gas-phase, UV-photo and ice 
chemistry, detailed heating & cooling balance, and continuum & 
line radiative transfer. 

The first results of the DENT grid show a strong dependence 
of the calculated emission line fluxes on the assumed stellar UV 
excess and on the flaring of the disc. The stellar UV is essential 
for the heating of the upper disc layers. In combination with posi- 
tive disc flaring, a strong stellar UV irradiation creates an extended 
warm surface layer with Tg > Td responsible for the line emissions. 
However, if the disc is not flared (self-shadowed), discs with total 
mass > 10"^ - 10"* Mo increasingly shield the stellar UV by 
their inner parts, which causes much cooler surface layers, and a 
saturation of the line fluxes with increasing disc mass. 

Despite these complicated parameter dependencies, we have 
shown that the [OI] 63.2/im line flux depends basically on two 
quantities, namely the total disc gas mass and the mean disc temper- 
ature. We will continue this work by two follow-up papers (Kamp 
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et al. 2010, in prep., Menard et al. 2010, in prep.) that will pro- 
vide more insight into the statistical behaviour of gas line and dust 
continuum predictions, respectively, to identify trends and robust 
correlations with disc mass. 

In summary, the DENT grid allows to 

• study the effects of stellar, disc, and dust parameters on contin- 
uum and Une observations, 

• allow for a quaUfied interpretation of observational data, 

• quickly predict line and continutmi fluxes for planning observa- 
tions, 

• search for best-fitting models concerning a given set of observed 
Une and continuum fluxes, 

• study the robustness of certain fit values against variation of the 
observational data. 

We intend to make the calculated DENT grid available to the scien- 
tific community. A graphical user interface called xDENT has been 
developed to allow researchers to visualise the DENT results, to 
make plots as presented in this letter, and to search for best-fitting 
models for a given set of continuum and line flux data. We em- 
phasise, however, that the DENT grid has not been developed for 
detailed fitting of individual objects. The coarse sampling of the 12- 
dimensional parameter space can mostly be used to narrow down 
the parameter range for individual objects, for example to design a 
finer sub-grid, especially for not so well-known objects. 

With a comprehensive data set of far IR gas emission lines to 
be obtained by Herschel/GASPS very soon, we aim at breaking the 
degeneracy of SED fitting and make possible a more profound anal- 
ysis of the physical, chemical and temperature structures of discs 
around young stars. 



APPENDIX A: SIMPLE LINE EMISSION MODEL 



Let us assume that an emission line is optically thin and that the 
emitting species is populated with a uniform excitation temperature 
The line luminosity is then given by 



iline = hi/ A^l Nu 



guexp(-_E„/fcTexc) 
' (3(rexc) 



(Al) 



where v is the line centre frequency, A^i is the Einstein coefficient 
of the line transition from level u to level I, Q is the partition func- 
tion, and Qu and Eu are the statistical weight and energy [K] of 
the upper level. The total number of line emitting particles iVtot is 
related to the total disc gas mass by 



(A2) 



where e is the abundance of the line emitting species with respect 
to hydrogen nuclei and fin 1.31 amu the gas mass per hydrogen 
nucleus, assuming solar abundances. 

The observable Une flux [erg/s/cm^] at distance d is 



-Fline — 



(A3) 



47rd2 ■ 

For the case of the [01] 63.2/im fine-structure line, we have u = 2, 
I = 1, A21 = 8.87 X 10-^ s-^, E2 = 227.7 K, 51 = 5 and 
g2 = 3. We calculate the partition function including the third level 
E3 — 326.0 K, = 1. The oxygen abundance assumed in the mod- 
els is e = 8.5 X 10"''. However, the actual abundance of the neutral 
oxygen atom is reduced by CO, CO-ice, H2O and H20-ice forma- 
tion, and we use a mean value from the models, £ = 3.6 x 10~*. 
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